1 Introducción

Los datos que se van a analizar en este proyecto han sido obtenidos desde Kaggle. Contienen precios de casas que fueron vendidas desde mayo de 2014 hasta mayo de 2015 en King County que es un condado ubicado en el estado estadounidense de Washington.

2 Objetivo del estudio

Lo que queremos hacer con estos datos es clasificar/predecir las viviendas conforme a su precio dependiendo de las variables recogidas de cada una. Para ello se implementarán varios modelos, comprobando cuáles de ellos funciona mejor con los tipos de variables que contamos. Finalmente se evaluarán los modelos en base a distintas métricas.

3 Datos

3.1 Categorización del precio

En nuestro estudio inicial, la variable “Precio” es continua, por lo que vamos a categorizarla. Para decidir las categorizaciones se ha usado los cuantiles. Se van a realizar dos tipos de categorizaciones:

  • Categorización 1: se ha categorizado en dos grupos, B1, casas baratas (casas con un precio < 500.000) y casas caras, C1, (casas con un precio > 500.000).
  • Categorización 2: se ha categorizado en tres grupos, B2, casas baratas (casas con un precio < 330.000), casas con un precio medio, M2 (330.000 < precio < 650.000) y casas caras, C2 (precio > 650.000).
#Categorizamos la variable respuesta price:
quantile(datos$price, prob=seq(0, 1, length = 5))
##      0%     25%     50%     75%    100% 
##   75000  321950  450000  645000 7700000
datos$price_categ1 <- cut(datos$price, breaks = c(0, 500000, 100000000), labels = c("B1", "C1"))
table(datos$price_categ1)
## 
##    B1    C1 
## 12560  9053
datos$price_categ2 <- cut(datos$price, breaks = c(0, 330000, 650000, 100000000), labels = c("B2","M2", "C2"))
table(datos$price_categ2)
## 
##    B2    M2    C2 
##  5881 10525  5207

En el siguiente mapa se puede visualizar cómo se distribuyen las casas “Caras” y “Baratas”. Se observa que las casas que están cercanas al agua y cerca de Seattle por la parte norte son más caras y hacia el sur son más baratas.

center_lon = median(datos$long,na.rm = TRUE)
center_lat = median(datos$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red"), 
                       datos$price_categ1 )

leaflet(datos) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal2(datos$price_categ1))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal2 , values = ~datos$price_categ1,
            title = "Tipos de Casas",
            opacity = 1)

En este otro mapa se puede visualizar cómo se distribuyen las casas según el precio en tres categorías:“Caras”, “Medio” y “Baratas”. Se puede observar cómo con esta categorización no está tan clara la separación entre clases.

center_lon = median(datos$long,na.rm = TRUE)
center_lat = median(datos$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red","yellow"), 
                       datos$price_categ2 )

leaflet(datos) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal2(datos$price_categ2))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal2 , values = ~datos$price_categ2,
            title = "Tipos de Casas 3 categorías",
            opacity = 1)

Por lo tanto, viendo ambos mapas, nos quedamos con la primera categorización: casas baratas y caras.

Eliminamos la segunda categorización:

datos$price_categ2= NULL

3.2 Train, test y validación

Se va a separar los datos en los 3 conjuntos de datos fundamentales:

  • Conjunto de datos de entrenamiento: en nuestro estudio datos_train, se corresponde con el 70% del total de los datos.
  • Conjunto de datos de validación: en nuestro estudio datos_validacion, se corresponde con el 15% del total de los datos.
  • Conjunto de datos de test: en nuestro estudio datos_test, se corresponde con el 15% del total de los datos.
num_total=nrow(datos)
set.seed(122556) #reproductividad

# 70% para train
indices_train = sample(1:num_total, .7*num_total)
datos_train = datos[indices_train,]

# 15% para test
indices=seq(1:num_total)
indices_test=indices[-indices_train]
indices_test1 = sample(indices_test, .15*num_total)
datos_test = datos[indices_test1,]

# 15% para validacion
indices_validacion=indices[c(-indices_train,-indices_test1)]
datos_validacion=datos[indices_validacion,]

3.3 Análisis exploratorio

Se van a realizar transformaciones de un conjunto de variables, estas transformaciones se aplicarán a cada conjunto de datos: train, test y validación.

  • Se realiza una transformación logarítmica sobre las variables sqft_living (pies cuadrados de la casa), sqft_lot (pies cuadrados del jardín) y sqft_above (pies cuadrados por encima del suelo). Hay que aclarar que esta última variable es la diferencia entre sqft_living y sqft_basement, por lo que va a estar altamentente correlada con sqft_living.

  • Se categorizan las variables:

    • Bathroom, esta varible puede tomar valores decimales de 0.25 en 0.25. El número de baños se contabiliza por las piezas y cada baño completo tiene 4 piezas. Por lo que con la nueva agrupación toma valores de 1 a 8 baños.

    • Sqft_basement, se categoriza como 0 las casas que no tienen sótano y 1 las casas que sí tienen sótano.

    • Grade, se va a categorizar del siguiente modo: con valor 0: calidad Baja, 1: calidad media y 2: calidad alta.

    • Year_renovated, se categoriza como 0: no ha tenido renovación y 1: sí ha tenido renovación.

  • Se pasan a factor las variables: waterfront, view, condition, grade_categ y zipcode.

  • Se eliminan Outliers.

3.3.1 Transformaciones datos Train

datos_train <- datos_train[,-2]

datos_train$id <- as.factor(datos_train$id)

datos_train$bathrooms_group <- cut(datos_train$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_train$bathrooms_group <- as.numeric(as.character(datos_train$bathrooms_group))

datos_train$log_sqft_living <- log10(datos_train$sqft_living)
datos_train$log_lot <- log10(datos_train$sqft_lot)
datos_train$log_above <- log10(datos_train$sqft_above)

datos_train$sqft_basement_cat <- cut(datos_train$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_train$waterfront<-as.factor(datos_train$waterfront)

datos_train$view<-as.factor(datos_train$view)

datos_train$condition<-as.factor(datos_train$condition)

datos_train$grade_categ <- cut(datos_train$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_train$yr_renovated_catg <-cut(datos_train$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_train$zipcode<-as.factor(datos_train$zipcode)

3.3.2 Eliminación de Outliers

datos_train$posicion<-c(1:nrow(datos_train))
indices_cero_habitaciones<-datos_train[datos_train$bedrooms==0,]$posicion
datos_train<-datos_train[-indices_cero_habitaciones,]

datos_train$posicion<-c(1:nrow(datos_train))
indices_cero_banos<-datos_train$posicion[datos_train$bathrooms_group==0]
datos_train<-datos_train[-indices_cero_banos,]

datos_train$posicion<-c(1:nrow(datos_train))
indice_hab33 <- datos_train[datos_train$bedrooms==33,]$posicion
datos_train[datos_train$posicion == indice_hab33,]$bedrooms = 3

3.3.3 Transformaciones adicionales:

Se han realizado dos categorizaciones adicionales sobre el conjunto de datos. Después de implementar varios modelos, se llegó a la conclusión de que algunas variables podían mejorar los resultados de los modelos siendo agrupadas. Para analizar cómo recategorizar estas variables se ha usado un árbol de decisión. Las variables son: zipcode y bathrooms_group.

  • Zipcode, esta variable es de tipo factor y tenía 70 códigos postales, por lo que se ha decidido aplicar un árbol de decisión para ver cómo clasificaba los códigos postales y así volver a categorizarla según el resultado obtenido.
model_selec_zipcode<-rpart(price_categ1~zipcode,data=datos_train ,parms=list(split="gini"))
print(model_selec_zipcode)
## n= 15119 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 15119 6385 B1 (0.5776837 0.4223163)  
##   2) zipcode=98001,98002,98003,98010,98011,98014,98019,98022,98023,98024,98028,98030,98031,98032,98034,98038,98042,98045,98055,98056,98058,98059,98070,98092,98106,98108,98118,98125,98126,98133,98144,98146,98148,98155,98166,98168,98178,98188,98198 8243 1336 B1 (0.8379231 0.1620769) *
##   3) zipcode=98004,98005,98006,98007,98008,98027,98029,98033,98039,98040,98052,98053,98065,98072,98074,98075,98077,98102,98103,98105,98107,98109,98112,98115,98116,98117,98119,98122,98136,98177,98199 6876 1827 C1 (0.2657068 0.7342932) *

Se va a categorizar en dos: Zona1 y Zona2.

datos_train$zona<-recode(datos_train$zipcode, "98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_train$zipcode = NULL
  • bathrooms_group, se aplica el mismo método que con zipcode para ver cómo se puede categorizar esta variable. Toma valores de 1 a 8 y queremos reducir el número de niveles.
model_selec_bathrooms<-rpart(price_categ1~bathrooms_group,data=datos_train ,parms=list(split="gini"))
print(model_selec_bathrooms)
## n= 15119 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 15119 6385 B1 (0.5776837 0.4223163)  
##   2) bathrooms_group< 2.5 7282 1850 B1 (0.7459489 0.2540511) *
##   3) bathrooms_group>=2.5 7837 3302 C1 (0.4213347 0.5786653) *

En el resultado del modelo se ve que corta en el número de baños < 2.5, por lo que se va a categorizar como 0 aquellas casas que tengan de 1 a 2 baños y como 1 las casas que tengan más de 2 baños.

datos_train$bathrooms_group <- cut(datos_train$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
# Limpiamos el dataframe
datos_train_limpio <- datos_train[c(3,21:25,8:10,26,16,17,29,27,20)]

#Eliminamos sqft_above porque es una combinalción lineal de sqft_living, están altamente correladas........

##### REVISAR COMENTARIO
datos_train_limpio$log_above = NULL

datos_train_numeric <- datos_train_limpio %>% select_if(is.numeric)

3.3.4 Transformaciones datos Test

Realizamos todas las transformaciones y categorizaciones para el conjunto de datos de test.

datos_test <- datos_test[,-2]

datos_test$id <- as.factor(datos_test$id)

datos_test$bathrooms_group <- cut(datos_test$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_test$bathrooms_group <- as.numeric(as.character(datos_test$bathrooms_group))

datos_test$log_sqft_living <- log10(datos_test$sqft_living)
datos_test$log_lot <- log10(datos_test$sqft_lot)
datos_test$log_above <- log10(datos_test$sqft_above)

datos_test$sqft_basement_cat <- cut(datos_test$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_test$waterfront<-as.factor(datos_test$waterfront)

datos_test$view<-as.factor(datos_test$view)

datos_test$condition<-as.factor(datos_test$condition)

datos_test$grade_categ <- cut(datos_test$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_test$yr_renovated_catg <-cut(datos_test$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_test$zipcode<-as.factor(datos_test$zipcode)

#codificar la variable Zipcode

datos_test$zona<-recode(datos_test$zipcode, " 98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_test$zipcode = NULL

datos_test$bathrooms_group <- cut(datos_test$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
datos_test_limpio <- datos_test[c(3,21:25,8:10,26,16,17,28,27,20)]
datos_test_limpio$log_above = NULL
datos_test_numeric <- datos_test_limpio %>% select_if(is.numeric)

3.3.5 Transformaciones datos Validación

Realizamos todas las transformaciones y categorizaciones para el conjunto de datos de Validación.

datos_validacion <- datos_validacion[,-2]

datos_validacion$id <- as.factor(datos_validacion$id)

datos_validacion$bathrooms_group <- cut(datos_validacion$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_validacion$bathrooms_group <- as.numeric(as.character(datos_validacion$bathrooms_group))

datos_validacion$log_sqft_living <- log10(datos_validacion$sqft_living)
datos_validacion$log_lot <- log10(datos_validacion$sqft_lot)
datos_validacion$log_above <- log10(datos_validacion$sqft_above)

datos_validacion$sqft_basement_cat <- cut(datos_validacion$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_validacion$waterfront<-as.factor(datos_validacion$waterfront)

datos_validacion$view<-as.factor(datos_validacion$view)

datos_validacion$condition<-as.factor(datos_validacion$condition)

datos_validacion$grade_categ <- cut(datos_validacion$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_validacion$yr_renovated_catg <-cut(datos_validacion$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_validacion$zipcode<-as.factor(datos_validacion$zipcode)

#codificar la variable Zipcode

datos_validacion$zona<-recode(datos_validacion$zipcode, " 98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_validacion$zipcode = NULL

datos_validacion$bathrooms_group <- cut(datos_validacion$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
datos_validacion_limpio <- datos_validacion[c(3,21:25,8:10,26,16,17,28,27,20)]
datos_validacion_limpio$log_above = NULL
datos_validacion_numeric <- datos_validacion_limpio %>% select_if(is.numeric)

3.3.6 Resumen EDA

La distribución de las casas según la nueva categorización daría como resultado un 57,8% de casas baratas (B1) y un 42,23% de casas caras (C1):

datos_train_limpio %>% 
  group_by(price_categ1) %>% 
  summarise(Count = n())%>% 
  mutate(percent = prop.table(Count)*100)%>%
  ggplot(aes(reorder(price_categ1, - percent), percent), fill = price_categ1)+
  geom_col(fill = c("red", "blue"))+
  xlab("Precio de las casas") + 
  ylab("Percent")+
  ggtitle("% de las casas por precio")

En cuanto a la relación de las variables que hacen referencia a las características de las casas con la variable respuesta categorizada, se observa lo siguiente:

p1<-ggplot(data = datos_train_limpio)+
  geom_bar(aes(x=bedrooms,fill=price_categ1,bins=30, position = "identity"))+
  facet_grid(price_categ1~., scales = 'free') 
p2<-ggplot(data=datos_train_limpio)+
  geom_bar(aes(x=bathrooms_group,fill = price_categ1))
p3<-ggplot(data=datos_train_limpio)+
  geom_density(aes(x=log_sqft_living, fill=price_categ1))
p4<-ggplot(data=datos_train_limpio)+
  geom_density(aes(x=log_lot, fill=price_categ1))+
  facet_grid(price_categ1~., scales = 'free')

grid.arrange(p1, p2, p3, p4, nrow = 2)

Por último, para la relación entre las variables que expresan la calidad de las casas y el precio, se obtienen los siguientes resultados:

p5<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=waterfront,fill=price_categ1, stat="count"))
p6<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=view,fill=price_categ1, stat="count"))
p7<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=condition,fill=price_categ1, stat="count"))
p8<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=grade_categ,fill=price_categ1, stat="count"))
p9<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=yr_renovated_catg,fill=price_categ1, stat="count"))
p10<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=zona,fill=price_categ1,bins=30, position = "identity"))

grid.arrange(p5,p6,p7,p8,p9,p10, nrow = 4)

4 Métodos de agrupamiento-No supervisado

A continuación, se implementan dos métodos de agrupamiento no supervisado. Primero, un método jerárquico para averiguar (si es posible) la cantidad adecuada de clústers y, a continuación, K-Means y K-Medoids.

4.1 Clúster Jerárquico

La agrupación es una técnica para agrupar puntos de datos similares en un grupo y separar las diferentes observaciones en diferentes grupos. En el Clustering Jerárquico los clústers se crean de manera que tengan un orden predeterminado. En nuestro estudio se va a aplicar un método aglomerativo que consiste en que cada observación se asigna a su propio clúster. Luego, se calcula la similitud (o distancia) entre cada uno de los clústers y los dos clústers más similares se fusionan en uno. Finalmente, los pasos 2 y 3 se repiten hasta que solo quede un grupo.

Aplicando este método a nuestros datos queremos observar si siguen algún patrón para poder agrupar las casas.

4.1.1 Train

Escalamos las variables numéricas, es decir, cada variable ahora tendrá una media cero y una desviación estándar uno. También se requieren los valores de distancia. La medida predeterminada para la función dist es ‘Euclidiana’.

datos_scale<- as.data.frame(scale(datos_train_numeric))

matriz_dist=dist(datos_scale)

En este caso particular, estimamos dos clústers:

modelo_hc= hclust(matriz_dist, method = "average")
plot(modelo_hc)
rect.hclust(modelo_hc, k=2,border="red")

A continuación vemos cómo se han agrupado los datos marcando que el número de clústers sea 2. Como se puede ver, prácticamente todas las casas están en el grupo 1.

grupos2=cutree(modelo_hc,k=2)
table(datos_train_limpio$price_categ1, grupos2)
##     grupos2
##         1    2
##   B1 8724   10
##   C1 6385    0

4.2 Clúster no Jerárquico- K-Means

El método de K-Means basa su funcionamiento en agrupar los datos de entrada en un total de k conjuntos definidos por un centroide, cuya distancia con los puntos que pertenecen a cada uno de los datos sea la menor posible.

4.2.1 Train

Primero se van a realizar los gráficos para ver cómo están diferenciadas las casas por precio. Para poder visualizarlo en dos dimensiones se ha usado la función “prcomp” (PCA)

colores1= c("red","blue")
colores11 = colores1[datos_train_limpio$price_categ1]

plot(prcomp(datos_train_numeric, scale = T)$x[,1:2], type="n",main= "Dos categorías")
text(prcomp(datos_train_numeric, scale = T)$x[,1:2], as.character(datos_train_limpio$price_categ1), col=colores11)

A continuación se va a aplicar el método de agrupamiento K-means, al igual que en el método anterior se le pasa la matriz de distacias y se van a agrupar los datos en dos conjuntos:

set.seed(1234)
model_km <- kmeans(matriz_dist, centers=2)
table(model_km$cluster) #asignación de observación a los cluster
## 
##     1     2 
##  3926 11193
colores1= c("red","blue")
colores11 = colores1[datos_train_limpio$price_categ1]

plot(prcomp(datos_train_numeric, scale = T)$x[,1:2], type="n")
text(prcomp(datos_train_numeric, scale = T)$x[,1:2], as.character(model_km$cluster), col=colores11)

4.2.1.1 Cálculo del k-óptimo.

Se va a determinar la cantidad óptima de centroides a utilizar a partir del Método del Codo. Para ello, aplicaremos la función kmeans al conjunto de datos, variando en cada caso el valor de k y acumulando los valores de WCSS (Within-Cluster-Sum-of-Squares) obtenidos:

set.seed(1234)
wcss <- vector()
for(i in 1:20){
  wcss[i] <- sum(kmeans(datos_scale, i)$withinss)
}

ggplot() + geom_point(aes(x = 1:20, y = wcss), color = 'blue') + 
  geom_line(aes(x = 1:20, y = wcss), color = 'blue') + 
  ggtitle("Método del Codo") + 
  xlab('Cantidad de Centroides k') + 
  ylab('WCSS')

Como se observa en la gráfica, el K-óptimo que se podría aplicar sería de 10.

A continuación, para visualizar si el agrupamiento que se ha llevado a cabo está relacionado con el precio de las casas (Caras, Baratas), se representa en el mapa que se muestra a continuación.

clusterkmeans=as.data.frame(model_km$cluster)
clusterkmeans$indice= as.integer(rownames(clusterkmeans))

colnames(clusterkmeans)[1]= "categoria_price_km"
clustering1= clusterkmeans[order(clusterkmeans$indice),]


center_lon = median(datos_train_limpio$long,na.rm = TRUE)
center_lat = median(datos_train_limpio$lat,na.rm = TRUE)

factpal1 <- colorFactor(c("green","red"),
                       clustering1$categoria_price_km )

leaflet(datos_train_limpio) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal1(clustering1$categoria_price_km))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal1 , values = ~clustering1$categoria_price_km,
            title = "Tipos de casas",
            opacity = 1)

Comparando el mapa con el de los datos iniciales, dista bastante. Por lo que deducimos que la agrupación que está realizando K-Means no es muy buena.

4.3 K-Medoids

K-medoids es un método de clustering muy similar a K-means en cuanto a que ambos agrupan las observaciones en K clusters, donde K es un valor preestablecido. La diferencia es que, en K-medoids, cada clúster está representado por una observación presente en el clúster (medoid). En nuestro estudio será una observación de una casa, mientras que en K-means cada clúster está representado por su centroide, que se corresponde con el promedio de todas las observaciones del clúster pero con ninguna casa en particular.

Este algoritmo es menos sensible al ruido y los valores atípicos, en comparación con k-means, porque usa medoides como centros de clúster en lugar de centroides. (utilizados en k-means).

datoskmedoids = datos_train_limpio[,-15]
model_medoids = pam(x = datoskmedoids, k = 2, keep.diss = TRUE, keep.data = TRUE)
model_medoids$medoids
##      bedrooms bathrooms_group log_sqft_living  log_lot sqft_basement_cat
## 9606        3               1        3.152288 3.870404                 1
## 8735        4               2        3.406540 3.936111                 2
##      waterfront view condition grade_categ     lat     long zona
## 9606          1    1         3           2 47.4949 -122.239    1
## 8735          1    1         3           2 47.6477 -122.197    2
##      yr_renovated_catg price_categ1
## 9606                 1            1
## 8735                 1            2
grupos<-data.frame(datoskmedoids)
grupos<-cbind(grupos,data.frame(model_medoids$clustering))
grupos$model_medoids.clustering<-as.factor(grupos$model_medoids.clustering)
ggpairs(grupos,columns=c(1,2,3,15,12),mapping=aes(color=model_medoids.clustering))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

table(grupos$model_medoids.clustering)
## 
##    1    2 
## 8831 6288
clustering= sort(model_medoids$clustering)
clustering=as.data.frame(model_medoids$clustering)
clustering$indice= as.integer(rownames(clustering))

colnames(clustering)[1]= "categoria_price"
clustering2= clustering[order(clustering$indice),]

center_lon = median(datoskmedoids$long,na.rm = TRUE)
center_lat = median(datoskmedoids$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red"), 
                       clustering2$categoria_price )

leaflet(datoskmedoids) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal2(clustering2$categoria_price))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal2 , values = ~clustering2$categoria_price,
            title = "Tipos de casas",
            opacity = 1)

5 Técnicas de reducción de la dimensionalidad

5.1 PCA

Es un método que permite simplificar la complejidad de espacios muestrales con muchas dimensiones a la vez que conserva su información.

pca<-prcomp(datos_train_numeric,scale=T)
plot(pca)

summary(pca)
## Importance of components:
##                          PC1    PC2    PC3    PC4     PC5
## Standard deviation     1.414 1.0885 0.9290 0.7850 0.57947
## Proportion of Variance 0.400 0.2369 0.1726 0.1232 0.06716
## Cumulative Proportion  0.400 0.6370 0.8096 0.9328 1.00000
biplot(x = pca, scale = 0, cex = 0.6, col = c("blue4", "brown3"))

6 Aprendizaje supervisado

6.1 GLM-Regresión Logística.

Con este modelo se va a estudiar si existe relación entre el hecho de que una casa sea “cara” ó “barata” dependiendo de las características de las casas. Se va a generar un modelo en el que a partir de las variables prediga la probabilidad de que una casa sea barata o cara.

6.1.1 Train

datos_train_rl <- datos_train_limpio[,-c(8,9)] # quitamos condition y grade_categ (no aportan **)
datos_train_rl$price_categ1<- recode(datos_train_rl$price_categ1, "'B1'=0; 'C1'=1")
model_glm = glm(price_categ1 ~., family = binomial, data =datos_train_rl )
summary(model_glm)
## 
## Call:
## glm(formula = price_categ1 ~ ., family = binomial, data = datos_train_rl)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.3536  -0.3998  -0.0883   0.3481   3.6096  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -759.91116   32.36457 -23.480  < 2e-16 ***
## bedrooms             -0.23062    0.03834  -6.015 1.80e-09 ***
## bathrooms_group1      0.14717    0.06828   2.155  0.03113 *  
## log_sqft_living      13.13156    0.32565  40.324  < 2e-16 ***
## log_lot               0.66619    0.08223   8.101 5.45e-16 ***
## sqft_basement_cat1   -0.45846    0.05867  -7.814 5.52e-15 ***
## waterfront1           1.35682    0.62801   2.161  0.03073 *  
## view1                 1.13634    0.22417   5.069 4.00e-07 ***
## view2                 1.21640    0.14137   8.605  < 2e-16 ***
## view3                 1.84563    0.21019   8.781  < 2e-16 ***
## view4                 2.49482    0.49846   5.005 5.58e-07 ***
## lat                   5.32574    0.24476  21.759  < 2e-16 ***
## long                 -3.75931    0.24366 -15.429  < 2e-16 ***
## zona2                 3.33880    0.07036  47.451  < 2e-16 ***
## yr_renovated_catg1    0.37667    0.13792   2.731  0.00631 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 20592.9  on 15118  degrees of freedom
## Residual deviance:  8975.8  on 15104  degrees of freedom
## AIC: 9005.8
## 
## Number of Fisher Scoring iterations: 6
# EVALUACION

predicciones <- ifelse(test = model_glm$fitted.values > 0.5, yes = 1, no = 0)
tabla_glm_train <- table(model_glm$model$price_categ1, predicciones,
                      dnn = c("observaciones", "predicciones"))
tabla_glm_train
##              predicciones
## observaciones    0    1
##             0 7813  921
##             1 1010 5375
#caras
pred_caras_glm_train=tabla_glm_train[2,2]/(tabla_glm_train[2,2]+tabla_glm_train[1,2])
rec_caras_glm_train=tabla_glm_train[2,2]/(tabla_glm_train[2,2]+tabla_glm_train[2,1])
F_caras_reglog=(2*pred_caras_glm_train*rec_caras_glm_train)/(pred_caras_glm_train+rec_caras_glm_train)
cat(c('F1 casas caras: ', F_caras_reglog), '\n')
## F1 casas caras:  0.847724942827853
#baratas
pred_baratas_glm_train=tabla_glm_train[1,1]/(tabla_glm_train[1,1]+tabla_glm_train[2,1])
rec_baratas_glm_train=tabla_glm_train[1,1]/(tabla_glm_train[1,1]+tabla_glm_train[1,2])
F_baratas_reglog=(2*pred_baratas_glm_train*rec_baratas_glm_train)/(pred_baratas_glm_train+rec_baratas_glm_train)
cat(c('F1 casas baratas: ', F_baratas_reglog), '\n')
## F1 casas baratas:  0.890015378481517
#F-Medida
F_reglog_train= (F_caras_reglog+F_baratas_reglog)/2
cat(c('F1 global: ', F_reglog_train), '\n')
## F1 global:  0.868870160654685
accuracy_reglog_train = (tabla_glm_train[1,1]+tabla_glm_train[2,2]) / (tabla_glm_train[1,1]+tabla_glm_train[1,2]+tabla_glm_train[2,1]+tabla_glm_train[2,2])
cat(c('Accuracy: ', accuracy_reglog_train), '\n')
## Accuracy:  0.872279912692638

6.2 KNN

6.2.1 Train

6.2.1.1 Ajuste de Hiperparámetros.

Hemos usado tres métodos para ajustar y comparar este modelo (parámetro k):

  1. Usando train.kknn.
  2. Usando tune.
  3. Usando diferentes k de 1 a 50 y comprobando los resultados de la F1-medida.

En primer lugar con la función train.kknn se obtiene de manera automática el mejor valor de k hasta un máximo de k=30.

set.seed(1234)

suppressWarnings(suppressMessages(library(kknn)))
knn_1 <- train.kknn(price_categ1 ~ ., data = datos_train_limpio, kmax = 30)
knn_1
## 
## Call:
## train.kknn(formula = price_categ1 ~ ., data = datos_train_limpio,     kmax = 30)
## 
## Type of response variable: nominal
## Minimal misclassification: 0.1149547
## Best kernel: optimal
## Best k: 15

En segundo lugar se usa la función tune que itera hasta k = 30 y muestra el error y la dispersión para cada valor de k.

set.seed(1234)

knn_2 <- tune.knn(x=scale(datos_train_numeric),
                  y=as.factor(datos_train_limpio$price_categ1), k = 1:30,
                  tunecontrol = tune.control(sampling = "boot"))
summary(knn_2)
## 
## Parameter tuning of 'knn.wrapper':
## 
## - sampling method: bootstrapping 
## 
## - best parameters:
##   k
##  27
## 
## - best performance: 0.1222577 
## 
## - Detailed performance results:
##     k     error  dispersion
## 1   1 0.1438725 0.003032594
## 2   2 0.1488328 0.004107803
## 3   3 0.1464719 0.003574503
## 4   4 0.1439387 0.004097428
## 5   5 0.1376046 0.003948778
## 6   6 0.1363546 0.003548451
## 7   7 0.1332349 0.003366625
## 8   8 0.1329151 0.003934370
## 9   9 0.1288270 0.002981856
## 10 10 0.1286652 0.003300452
## 11 11 0.1273676 0.003302001
## 12 12 0.1269916 0.004677623
## 13 13 0.1258752 0.004332600
## 14 14 0.1269473 0.003960220
## 15 15 0.1262791 0.003599693
## 16 16 0.1249677 0.004310555
## 17 17 0.1255292 0.004179737
## 18 18 0.1256711 0.003900466
## 19 19 0.1251843 0.003735257
## 20 20 0.1245178 0.003857497
## 21 21 0.1247299 0.003625392
## 22 22 0.1252001 0.003437002
## 23 23 0.1242136 0.003550403
## 24 24 0.1235049 0.003724575
## 25 25 0.1230196 0.004134106
## 26 26 0.1224334 0.003672093
## 27 27 0.1222577 0.004162912
## 28 28 0.1224645 0.003954248
## 29 29 0.1229364 0.003490636
## 30 30 0.1235549 0.002958481
plot(knn_2, main="Mejor k según tune")

En tercer lugar, se comprueba manualmente los resultados de la F1-medida con diferentes modelos desde \(k=1\) hasta el \(k=50\) y se representa el resultado.

set.seed(1234)

k_maximo=50

rango=1:k_maximo
f1_modelos=c()

for (i in rango){
  model_knn=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=i)
  tabla=table(datos_train_limpio$price_categ1,model_knn)
  # f1 casas caras
  pred_means_caras=tabla[2,2]/(tabla[2,2]+tabla[1,2])
  rec_means_caras=tabla[2,2]/(tabla[2,2]+tabla[2,1])
  f1_caras=(2*pred_means_caras*rec_means_caras)/(pred_means_caras+rec_means_caras)
  # f1 casas baratas
  pred_means_baratas=tabla[1,1]/(tabla[1,1]+tabla[2,1])
  rec_means_baratas=tabla[1,1]/(tabla[1,1]+tabla[1,2])
  f1_baratas=(2*pred_means_baratas*rec_means_baratas)/(pred_means_baratas+rec_means_baratas)
  f1_total = (f1_baratas + f1_caras)/2
  f1_modelos=c(f1_modelos,f1_total)
}


plot(f1_modelos)

cat(c('Valor óptimo de k: ', which.max(f1_modelos)))
## Valor óptimo de k:  17

De la primer forma obtenemos un valor óptimo de \(k=15\), del segundo modo \(k=27\), y de la forma manual, concluimos que el mejor valor de k en función de la F1-medida es de \(k=17\). Finalmente, elegimos el valor de \(k=17\), implementamos el modelo y lo evaluamos.

model_knn=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=17, prob = TRUE)
center_lon = median(datos_train_limpio$long,na.rm = TRUE)
center_lat = median(datos_train_limpio$lat,na.rm = TRUE)

factpal1 <- colorFactor(c("green","red"), 
                       model_knn)

leaflet(datos_train_limpio) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal1(model_knn))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 12) %>%

  addLegend("bottomright", pal = factpal1 , values = ~model_knn,
            title = "Precio (en miles de $)",
            opacity = 1)
tabla_knn_train=table(datos_train_limpio$price_categ1, model_knn)
tabla_knn_train
##     model_knn
##        B1   C1
##   B1 7876  858
##   C1  915 5470
pred_caras_knn_train = tabla_knn_train[2,2]/(tabla_knn_train[2,2]+tabla_knn_train[1,2])
rec_caras_knn_train = tabla_knn_train[2,2]/(tabla_knn_train[2,2]+tabla_knn_train[2,1])
F_caras_knn_train=(2*pred_caras_knn_train*rec_caras_knn_train)/(pred_caras_knn_train+rec_caras_knn_train)
cat(c('F1 caras: ', F_caras_knn_train), '\n')
## F1 caras:  0.860536458743019
pred_baratas_knn_train=tabla_knn_train[1,1]/(tabla_knn_train[1,1]+tabla_knn_train[2,1])
rec_baratas_knn_train=tabla_knn_train[1,1]/(tabla_knn_train[1,1]+tabla_knn_train[1,2])
F_baratas_knn_train=(2*pred_baratas_knn_train*rec_baratas_knn_train)/(pred_baratas_knn_train+rec_baratas_knn_train)
cat(c('F1 baratas: ', F_baratas_knn_train), '\n')
## F1 baratas:  0.898830242510699
F_knn_train= (F_caras_knn_train+F_baratas_knn_train)/2
cat(c('F1 global: ', F_knn_train), '\n')
## F1 global:  0.879683350626859
accuracy_knn_train= (tabla_knn_train[1,1]+tabla_knn_train[2,2]) / (tabla_knn_train[1,1]+tabla_knn_train[1,2]+tabla_knn_train[2,1]+tabla_knn_train[2,2])
cat(c('Accuracy: ', accuracy_knn_train), '\n')
## Accuracy:  0.882730339308155

6.2.2 Test

tabla_knn_test=table(datos_test_limpio$price_categ1, knn(scale(datos_train_numeric), scale(datos_test_numeric),cl=datos_train_limpio$price_categ1,k=17))
tabla_knn_test
##     
##        B1   C1
##   B1 1701  210
##   C1  187 1143
pred_caras_knn_test=tabla_knn_test[2,2]/(tabla_knn_test[2,2]+tabla_knn_test[1,2])
rec_caras_knn_test=tabla_knn_test[2,2]/(tabla_knn_test[2,2]+tabla_knn_test[2,1])
F_caras_knn_test=(2*pred_caras_knn_test*rec_caras_knn_test)/(pred_caras_knn_test+rec_caras_knn_test)
cat(c('F1 caras: ', F_caras_knn_test), '\n')
## F1 caras:  0.852031308237048
pred_baratas_knn_test=tabla_knn_test[1,1]/(tabla_knn_test[1,1]+tabla_knn_test[2,1])
rec_baratas_knn_test=tabla_knn_test[1,1]/(tabla_knn_test[1,1]+tabla_knn_test[1,2])
F_baratas_knn_test=(2*pred_baratas_knn_test*rec_baratas_knn_test)/(pred_baratas_knn_test+rec_baratas_knn_test)
cat(c('F1 baratas: ', F_baratas_knn_test), '\n')
## F1 baratas:  0.895498815477757
F_knn_test = (F_caras_knn_test+F_baratas_knn_test)/2
cat(c('F1 global: ', F_knn_test), '\n')
## F1 global:  0.873765061857403
accuracy_knn_test= (tabla_knn_test[1,1]+tabla_knn_test[2,2]) / (tabla_knn_test[1,1]+tabla_knn_test[1,2]+tabla_knn_test[2,1]+tabla_knn_test[2,2])
cat(c('Accuracy: ', accuracy_knn_test), '\n')
## Accuracy:  0.877506942301759

6.2.3 Validación

tabla_knn_validacion=table(datos_validacion_limpio$price_categ1, knn(scale(datos_train_numeric), scale(datos_validacion_numeric),cl=datos_train_limpio$price_categ1,k=17))
tabla_knn_validacion
##     
##        B1   C1
##   B1 1727  179
##   C1  203 1134
pred_caras_knn_validacion=tabla_knn_validacion[2,2]/(tabla_knn_validacion[2,2]+tabla_knn_validacion[1,2])
rec_caras_knn_validacion=tabla_knn_validacion[2,2]/(tabla_knn_validacion[2,2]+tabla_knn_validacion[2,1])
F_caras_knn_validacion=(2*pred_caras_knn_validacion*rec_caras_knn_validacion)/(pred_caras_knn_validacion+rec_caras_knn_validacion)
cat(c('F1 caras: ', F_caras_knn_validacion), '\n')
## F1 caras:  0.855849056603774
pred_baratas_knn_validacion=tabla_knn_validacion[1,1]/(tabla_knn_validacion[1,1]+tabla_knn_validacion[2,1])
rec_baratas_knn_validacion=tabla_knn_validacion[1,1]/(tabla_knn_validacion[1,1]+tabla_knn_validacion[1,2])
F_baratas_knn_validacion=(2*pred_baratas_knn_validacion*rec_baratas_knn_validacion)/(pred_baratas_knn_validacion+rec_baratas_knn_validacion)
cat(c('F1 baratas: ', F_baratas_knn_validacion), '\n')
## F1 baratas:  0.900417101147028
F_knn_validacion = (F_caras_knn_validacion+F_baratas_knn_validacion)/2
cat(c('F1 global: ', F_knn_validacion), '\n')
## F1 global:  0.878133078875401
accuracy_knn_validacion= (tabla_knn_validacion[1,1]+tabla_knn_validacion[2,2]) / (tabla_knn_validacion[1,1]+tabla_knn_validacion[1,2]+tabla_knn_validacion[2,1]+tabla_knn_validacion[2,2])
cat(c('Accuracy: ', accuracy_knn_validacion), '\n')
## Accuracy:  0.882207832254086

6.3 Decision Trees

6.3.1 Train

#quitamos lat y long:
datos_decision_tree<-datos_train_limpio[,-c(10,11)]


decisiontree_model=rpart(price_categ1~., 
                         data=datos_decision_tree, 
                         parms=list(split="gini"),
                         control = rpart.control(cp = 0.005, maxdepth = 5, minbucket = 400))
# print(decisiontree_model)

fancyRpartPlot(decisiontree_model, caption='')

tabla_train_arbol=table(obs = datos_decision_tree$price_categ1, pred = predict(decisiontree_model, type = "class"))
tabla_train_arbol
##     pred
## obs    B1   C1
##   B1 7896  838
##   C1 1461 4924
pred_caras_dt_train = tabla_train_arbol[2,2]/(tabla_train_arbol[2,2]+tabla_train_arbol[1,2])
rec_caras_dt_train = tabla_train_arbol[2,2]/(tabla_train_arbol[2,2]+tabla_train_arbol[2,1])
F_caras_dt_train=(2*pred_caras_dt_train*rec_caras_dt_train)/(pred_caras_dt_train+rec_caras_dt_train)
cat(c('F1 caras: ', F_caras_dt_train), '\n')
## F1 caras:  0.810735160945089
pred_baratas_dt_train = tabla_train_arbol[1,1]/(tabla_train_arbol[1,1]+tabla_train_arbol[2,1])
rec_baratas_dt_train = tabla_train_arbol[1,1]/(tabla_train_arbol[1,1]+tabla_train_arbol[1,2])
F_baratas_dt_train=(2*pred_baratas_dt_train*rec_baratas_dt_train)/(pred_baratas_dt_train+rec_baratas_dt_train)
cat(c('F1 baratas: ', F_baratas_dt_train), '\n')
## F1 baratas:  0.872920236581726
F_dt_train= (F_caras_dt_train+F_baratas_dt_train)/2
cat(c('F1 global: ', F_dt_train), '\n')
## F1 global:  0.841827698763407
accuracy_dt_train = (tabla_train_arbol[1,1]+tabla_train_arbol[2,2]) / (tabla_train_arbol[1,1]+tabla_train_arbol[1,2]+tabla_train_arbol[2,1]+tabla_train_arbol[2,2])
cat(c('Accuracy: ', F_dt_train), '\n')
## Accuracy:  0.841827698763407

6.3.2 Test

tabla_test_arbol=table(obs = datos_test_limpio$price_categ1, pred = predict(decisiontree_model, datos_test_limpio[,-15], type = "class"))
tabla_test_arbol
##     pred
## obs    B1   C1
##   B1 1699  212
##   C1  302 1028
pred_caras_dt_test = tabla_test_arbol[2,2]/(tabla_test_arbol[2,2]+tabla_test_arbol[1,2])
rec_caras_dt_test = tabla_test_arbol[2,2]/(tabla_test_arbol[2,2]+tabla_test_arbol[2,1])
F_caras_dt_test=(2*pred_caras_dt_test*rec_caras_dt_test)/(pred_caras_dt_test+rec_caras_dt_test)
cat(c('F1 caras: ', F_caras_dt_test), '\n')
## F1 caras:  0.8
pred_baratas_dt_test = tabla_test_arbol[1,1]/(tabla_test_arbol[1,1]+tabla_test_arbol[2,1])
rec_baratas_dt_test = tabla_test_arbol[1,1]/(tabla_test_arbol[1,1]+tabla_test_arbol[1,2])
F_baratas_dt_test=(2*pred_baratas_dt_test*rec_baratas_dt_test)/(pred_baratas_dt_test+rec_baratas_dt_test)
cat(c('F1 baratas: ', F_baratas_dt_test), '\n')
## F1 baratas:  0.868609406952965
F_dt_test= (F_caras_dt_test+F_baratas_dt_test)/2
cat(c('F1 global: ', F_dt_test), '\n')
## F1 global:  0.834304703476483
accuracy_dt_test = (tabla_test_arbol[1,1]+tabla_test_arbol[2,2]) / (tabla_test_arbol[1,1]+tabla_test_arbol[1,2]+tabla_test_arbol[2,1]+tabla_test_arbol[2,2])
cat(c('Accuracy: ', F_dt_test), '\n')
## Accuracy:  0.834304703476483

6.3.3 Validación

tabla_validacion_arbol=table(obs = datos_validacion_limpio$price_categ1, pred = predict(decisiontree_model, datos_validacion_limpio[,-15], type = "class"))
tabla_validacion_arbol
##     pred
## obs    B1   C1
##   B1 1709  197
##   C1  301 1036
pred_caras_dt_validacion = tabla_validacion_arbol[2,2]/(tabla_validacion_arbol[2,2]+tabla_validacion_arbol[1,2])
rec_caras_dt_validacion = tabla_validacion_arbol[2,2]/(tabla_validacion_arbol[2,2]+tabla_validacion_arbol[2,1])
F_caras_dt_validacion=(2*pred_caras_dt_validacion*rec_caras_dt_validacion)/(pred_caras_dt_validacion+rec_caras_dt_validacion)
cat(c('F1 caras: ', F_caras_dt_validacion), '\n')
## F1 caras:  0.806225680933852
pred_baratas_dt_validacion = tabla_validacion_arbol[1,1]/(tabla_validacion_arbol[1,1]+tabla_validacion_arbol[2,1])
rec_baratas_dt_validacion = tabla_validacion_arbol[1,1]/(tabla_validacion_arbol[1,1]+tabla_validacion_arbol[1,2])
F_baratas_dt_validacion=(2*pred_baratas_dt_validacion*rec_baratas_dt_validacion)/(pred_baratas_dt_validacion+rec_baratas_dt_validacion)
cat(c('F1 baratas: ', F_baratas_dt_validacion), '\n')
## F1 baratas:  0.872829417773238
F_dt_validacion= (F_caras_dt_validacion+F_baratas_dt_validacion)/2
cat(c('F1 global: ', F_dt_validacion), '\n')
## F1 global:  0.839527549353545
accuracy_dt_validacion = (tabla_validacion_arbol[1,1]+tabla_validacion_arbol[2,2]) / (tabla_validacion_arbol[1,1]+tabla_validacion_arbol[1,2]+tabla_validacion_arbol[2,1]+tabla_validacion_arbol[2,2])
cat(c('Accuracy: ', F_dt_validacion), '\n')
## Accuracy:  0.839527549353545

6.4 Random Forest

6.4.1 Train

set.seed(1234)
randomforest_model=randomForest(price_categ1~., 
                                data=datos_train_limpio, 
                                parms=list(split="gini"),
                                ntree=200, 
                                importance = FALSE, 
                                proximity = FALSE, 
                                mtry=6, 
                                replace = TRUE)
print(randomforest_model)
## 
## Call:
##  randomForest(formula = price_categ1 ~ ., data = datos_train_limpio,      parms = list(split = "gini"), ntree = 200, importance = FALSE,      proximity = FALSE, mtry = 6, replace = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 200
## No. of variables tried at each split: 6
## 
##         OOB estimate of  error rate: 9.56%
## Confusion matrix:
##      B1   C1 class.error
## B1 7999  735  0.08415388
## C1  710 5675  0.11119812
tabla_train_randomforest = table(obs = datos_train_limpio$price_categ1, pred = predict(randomforest_model, type = "class") )
tabla_train_randomforest
##     pred
## obs    B1   C1
##   B1 7999  735
##   C1  710 5675
pred_caras_rf_train = tabla_train_randomforest[2,2]/(tabla_train_randomforest[2,2]+tabla_train_randomforest[1,2])
rec_caras_rf_train = tabla_train_randomforest[2,2]/(tabla_train_randomforest[2,2]+tabla_train_randomforest[2,1])
F_caras_rf_train=(2*pred_caras_rf_train*rec_caras_rf_train)/(pred_caras_rf_train+rec_caras_rf_train)
cat(c('F1 caras: ', F_caras_rf_train), '\n')
## F1 caras:  0.887065259867136
pred_baratas_rf_train = tabla_train_randomforest[1,1]/(tabla_train_randomforest[1,1]+tabla_train_randomforest[2,1])
rec_baratas_rf_train=tabla_train_randomforest[1,1]/(tabla_train_randomforest[1,1]+tabla_train_randomforest[1,2])
F_baratas_rf_train=(2*pred_baratas_rf_train*rec_baratas_rf_train)/(pred_baratas_rf_train+rec_baratas_rf_train)
cat(c('F1 baratas: ', F_baratas_rf_train), '\n')
## F1 baratas:  0.917158745628619
F_rf_train= (F_caras_rf_train+F_baratas_rf_train)/2
cat(c('F1 global: ', F_rf_train), '\n')
## F1 global:  0.902112002747877
accuracy_rf_train = (tabla_train_randomforest[1,1]+tabla_train_randomforest[2,2]) / (tabla_train_randomforest[1,1]+tabla_train_randomforest[1,2]+tabla_train_randomforest[2,1]+tabla_train_randomforest[2,2])
cat(c('Accuracy: ', accuracy_rf_train), '\n')
## Accuracy:  0.904424895826444

6.4.2 Test

tabla_test_randomforest = table(obs = datos_test_limpio$price_categ1, pred = predict(randomforest_model, datos_test_limpio[,-15], type = "class") )
tabla_test_randomforest
##     pred
## obs    B1   C1
##   B1 1751  160
##   C1  138 1192
pred_caras_rf_test = tabla_test_randomforest[2,2]/(tabla_test_randomforest[2,2]+tabla_test_randomforest[1,2])
rec_caras_rf_test = tabla_test_randomforest[2,2]/(tabla_test_randomforest[2,2]+tabla_test_randomforest[2,1])
F_caras_rf_test=(2*pred_caras_rf_test*rec_caras_rf_test)/(pred_caras_rf_test+rec_caras_rf_test)
cat(c('F1 caras: ', F_caras_rf_test), '\n')
## F1 caras:  0.888888888888889
pred_baratas_rf_test = tabla_test_randomforest[1,1]/(tabla_test_randomforest[1,1]+tabla_test_randomforest[2,1])
rec_baratas_rf_test=tabla_test_randomforest[1,1]/(tabla_test_randomforest[1,1]+tabla_test_randomforest[1,2])
F_baratas_rf_test=(2*pred_baratas_rf_test*rec_baratas_rf_test)/(pred_baratas_rf_test+rec_baratas_rf_test)
cat(c('F1 baratas: ', F_baratas_rf_test), '\n')
## F1 baratas:  0.921578947368421
F_rf_test= (F_caras_rf_test+F_baratas_rf_test)/2
cat(c('F1 global: ', F_rf_test), '\n')
## F1 global:  0.905233918128655
accuracy_rf_test = (tabla_test_randomforest[1,1]+tabla_test_randomforest[2,2]) / (tabla_test_randomforest[1,1]+tabla_test_randomforest[1,2]+tabla_test_randomforest[2,1]+tabla_test_randomforest[2,2])
cat(c('Accuracy: ', accuracy_rf_test), '\n')
## Accuracy:  0.908053070040111

6.4.3 Validación

tabla_validacion_randomforest = table(obs = datos_validacion_limpio$price_categ1, pred = predict(randomforest_model, datos_validacion_limpio[,-15], type = "class") )
tabla_validacion_randomforest
##     pred
## obs    B1   C1
##   B1 1755  151
##   C1  147 1190
pred_caras_rf_validacion = tabla_validacion_randomforest[2,2]/(tabla_validacion_randomforest[2,2]+tabla_validacion_randomforest[1,2])
rec_caras_rf_validacion = tabla_validacion_randomforest[2,2]/(tabla_validacion_randomforest[2,2]+tabla_validacion_randomforest[2,1])
F_caras_rf_validacion=(2*pred_caras_rf_validacion*rec_caras_rf_validacion)/(pred_caras_rf_validacion+rec_caras_rf_validacion)
cat(c('F1 caras: ', F_caras_rf_validacion), '\n')
## F1 caras:  0.888722927557879
pred_baratas_rf_validacion = tabla_validacion_randomforest[1,1]/(tabla_validacion_randomforest[1,1]+tabla_validacion_randomforest[2,1])
rec_baratas_rf_validacion=tabla_validacion_randomforest[1,1]/(tabla_validacion_randomforest[1,1]+tabla_validacion_randomforest[1,2])
F_baratas_rf_validacion=(2*pred_baratas_rf_validacion*rec_baratas_rf_validacion)/(pred_baratas_rf_validacion+rec_baratas_rf_validacion)
cat(c('F1 baratas: ', F_baratas_rf_validacion), '\n')
## F1 baratas:  0.921743697478992
F_rf_validacion= (F_caras_rf_validacion+F_baratas_rf_validacion)/2
cat(c('F1 global: ', F_rf_validacion), '\n')
## F1 global:  0.905233312518435
accuracy_rf_validacion = (tabla_validacion_randomforest[1,1]+tabla_validacion_randomforest[2,2]) / (tabla_validacion_randomforest[1,1]+tabla_validacion_randomforest[1,2]+tabla_validacion_randomforest[2,1]+tabla_validacion_randomforest[2,2])
cat(c('Accuracy: ', accuracy_rf_validacion), '\n')
## Accuracy:  0.908109774899784

6.5 SVM

set.seed(1234)
modelo_svm <- e1071::svm(formula = price_categ1 ~., 
                         data = datos_train_limpio, 
                         kernel = "linear", 
                         probability =TRUE,
                         cost=1)
summary(modelo_svm)
## 
## Call:
## svm(formula = price_categ1 ~ ., data = datos_train_limpio, kernel = "linear", 
##     probability = TRUE, cost = 1)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1 
## 
## Number of Support Vectors:  4702
## 
##  ( 2352 2350 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  B1 C1

6.5.1 Train

tabla_svm_train=table(obs = datos_train_limpio$price_categ1, pred = predict(modelo_svm,datos_train_limpio[,-15], type ="class"))
tabla_svm_train
##     pred
## obs    B1   C1
##   B1 7835  899
##   C1 1034 5351
pred_caras_svm_train=tabla_svm_train[2,2]/(tabla_svm_train[2,2]+tabla_svm_train[1,2])
rec_caras_svm_train=tabla_svm_train[2,2]/(tabla_svm_train[2,2]+tabla_svm_train[2,1])
F_caras_svm_train=(2*pred_caras_svm_train*rec_caras_svm_train)/(pred_caras_svm_train+rec_caras_svm_train)
cat(c('F1 caras: ', F_caras_svm_train), '\n')
## F1 caras:  0.847012267510883
pred_baratas_svm_train=tabla_svm_train[1,1]/(tabla_svm_train[1,1]+tabla_svm_train[2,1])
rec_baratas_svm_train=tabla_svm_train[1,1]/(tabla_svm_train[1,1]+tabla_svm_train[1,2])
F_baratas_svm_train=(2*pred_baratas_svm_train*rec_baratas_svm_train)/(pred_baratas_svm_train+rec_baratas_svm_train)
cat(c('F1 baratas: ', F_baratas_svm_train), '\n')
## F1 baratas:  0.890189172300176
F_svm_train= (F_caras_svm_train+F_baratas_svm_train)/2
cat(c('F1 global: ', F_svm_train), '\n')
## F1 global:  0.868600719905529
accuracy_svm_train = (tabla_svm_train[1,1]+tabla_svm_train[2,2]) / (tabla_svm_train[1,1]+tabla_svm_train[1,2]+tabla_svm_train[2,1]+tabla_svm_train[2,2])
cat(c('Accuracy: ', accuracy_svm_train), '\n')
## Accuracy:  0.872147628811429

6.5.2 Test

tabla_svm_test=table(obs = datos_test_limpio$price_categ1,  pred = predict(modelo_svm,datos_test_limpio[,-15], type ="class"))
tabla_svm_test
##     pred
## obs    B1   C1
##   B1 1694  217
##   C1  213 1117
pred_caras_svm_test=tabla_svm_test[2,2]/(tabla_svm_test[2,2]+tabla_svm_test[1,2])
rec_caras_svm_test=tabla_svm_test[2,2]/(tabla_svm_test[2,2]+tabla_svm_test[2,1])
F_caras_svm_test=(2*pred_caras_svm_test*rec_caras_svm_test)/(pred_caras_svm_test+rec_caras_svm_test)
cat(c('F1 caras: ', F_caras_svm_test), '\n')
## F1 caras:  0.838588588588589
pred_baratas_svm_test=tabla_svm_test[1,1]/(tabla_svm_test[1,1]+tabla_svm_test[2,1])
rec_baratas_svm_test=tabla_svm_test[1,1]/(tabla_svm_test[1,1]+tabla_svm_test[1,2])
F_baratas_svm_test=(2*pred_baratas_svm_test*rec_baratas_svm_test)/(pred_baratas_svm_test+rec_baratas_svm_test)
cat(c('F1 baratas: ', F_baratas_svm_test), '\n')
## F1 baratas:  0.887375589313777
F_svm_test= (F_caras_svm_test+F_baratas_svm_test)/2
cat(c('F1 global: ', F_svm_test), '\n')
## F1 global:  0.862982088951183
accuracy_svm_test = (tabla_svm_test[1,1]+tabla_svm_test[2,2]) / (tabla_svm_test[1,1]+tabla_svm_test[1,2]+tabla_svm_test[2,1]+tabla_svm_test[2,2])
cat(c('Accuracy: ', accuracy_svm_test), '\n')
## Accuracy:  0.867324899722308

6.5.3 Validación

tabla_svm_validacion=table(obs = datos_validacion_limpio$price_categ1,  pred = predict(modelo_svm,datos_validacion_limpio[,-15], type ="class"))
tabla_svm_validacion
##     pred
## obs    B1   C1
##   B1 1711  195
##   C1  206 1131
pred_caras_svm_validacion=tabla_svm_validacion[2,2]/(tabla_svm_validacion[2,2]+tabla_svm_validacion[1,2])
rec_caras_svm_validacion=tabla_svm_validacion[2,2]/(tabla_svm_validacion[2,2]+tabla_svm_validacion[2,1])
F_caras_svm_validacion=(2*pred_caras_svm_validacion*rec_caras_svm_validacion)/(pred_caras_svm_validacion+rec_caras_svm_validacion)
cat(c('F1 caras: ', F_caras_svm_validacion), '\n')
## F1 caras:  0.849417949680811
pred_baratas_svm_validacion=tabla_svm_validacion[1,1]/(tabla_svm_validacion[1,1]+tabla_svm_validacion[2,1])
rec_baratas_svm_validacion=tabla_svm_validacion[1,1]/(tabla_svm_validacion[1,1]+tabla_svm_validacion[1,2])
F_baratas_svm_validacion=(2*pred_baratas_svm_validacion*rec_baratas_svm_validacion)/(pred_baratas_svm_validacion+rec_baratas_svm_validacion)
cat(c('F1 baratas: ', F_baratas_svm_validacion), '\n')
## F1 baratas:  0.895108553492022
F_svm_validacion= (F_caras_svm_validacion+F_baratas_svm_validacion)/2
cat(c('F1 global: ', F_svm_validacion), '\n')
## F1 global:  0.872263251586417
accuracy_svm_validacion = (tabla_svm_validacion[1,1]+tabla_svm_validacion[2,2]) / (tabla_svm_validacion[1,1]+tabla_svm_validacion[1,2]+tabla_svm_validacion[2,1]+tabla_svm_validacion[2,2])
cat(c('Accuracy: ', accuracy_svm_validacion), '\n')
## Accuracy:  0.876349059512797

6.5.4 Ajuste de hiperparámetros

svm_tune <- tune(svm, price_categ1 ~., data = datos_train_limpio, ranges = list(cost = 2^(2:4)), tunecontrol = tune.control(sampling = "fix"))
summary(svm_tune)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: fixed training/validation set 
## 
## - best parameters:
##  cost
##    16
## 
## - best performance: 0.110119 
## 
## - Detailed performance results:
##   cost     error dispersion
## 1    4 0.1142857         NA
## 2    8 0.1125000         NA
## 3   16 0.1101190         NA
plot(svm_tune)

Pese a que el ajuste de hiperparámetros indica que el valor óptimo de cost es 8, repetimos los cálculos y la F1-global mejora muy poco (aprox. 0.0004). Pasamos de una F1 de 0.8686 a una F1 de 0.8690.

7 GAM

7.0.1 Train

datos_train_gam <- datos_train_limpio
datos_train_gam$price <- datos_train$price
 
model_gam <- gam(price ~ s(bedrooms, bs = "cr") + bathrooms_group + s(log_sqft_living, bs = "cr") + 
                   s(log_lot, bs = "cr") + sqft_basement_cat + waterfront + view + s(lat, bs = "cr") + 
                   s(long, bs = "cr") + zona + yr_renovated_catg, data = datos_train_gam)

summary(model_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## price ~ s(bedrooms, bs = "cr") + bathrooms_group + s(log_sqft_living, 
##     bs = "cr") + s(log_lot, bs = "cr") + sqft_basement_cat + 
##     waterfront + view + s(lat, bs = "cr") + s(long, bs = "cr") + 
##     zona + yr_renovated_catg
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          460094       3772 121.978  < 2e-16 ***
## bathrooms_group1      35433       4091   8.662  < 2e-16 ***
## sqft_basement_cat1   -33249       3315 -10.029  < 2e-16 ***
## waterfront1          523178      20514  25.503  < 2e-16 ***
## view1                 90905      11827   7.686 1.61e-14 ***
## view2                 82058       7077  11.595  < 2e-16 ***
## view3                155689       9712  16.031  < 2e-16 ***
## view4                287662      14672  19.606  < 2e-16 ***
## zona2                126977       5000  25.398  < 2e-16 ***
## yr_renovated_catg1    58719       7070   8.306  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df      F p-value    
## s(bedrooms)        8.739  8.957   19.8  <2e-16 ***
## s(log_sqft_living) 8.911  8.997 1509.6  <2e-16 ***
## s(log_lot)         8.339  8.802   37.3  <2e-16 ***
## s(lat)             8.872  8.995  280.6  <2e-16 ***
## s(long)            8.935  8.999  171.2  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.781   Deviance explained = 78.2%
## GCV = 3.0207e+10  Scale est. = 3.0099e+10  n = 15119
plot(model_gam, ylab = "price")

par(mfrow=c(3,3))
visreg(model_gam)

7.0.2 Test

datos_test_limpio$price <- datos_test$price
datos_test_limpio$price_predict = predict(model_gam, datos_test_limpio)


datos_test_limpio$precio_diff <- abs(datos_test_limpio$price - datos_test_limpio$price_predict)

dif_media_test <- mean(datos_test_limpio$precio_diff)
dif_media_test
## [1] 108157.9

7.0.3 Validación

datos_validacion_limpio$price <- datos_validacion$price
datos_validacion_limpio$price_predict = predict(model_gam, datos_validacion_limpio)


datos_validacion_limpio$precio_diff <- abs(datos_validacion_limpio$price - datos_validacion_limpio$price_predict)

dif_media_test <- mean(datos_validacion_limpio$precio_diff)
dif_media_test
## [1] 105924.8

8 Evaluación y comparación de modelos

Una vez que se han entrenado y optimizado distintos modelos, se tiene que identificar cuál de ellos consigue mejores resultados para el problema en cuestión, en este caso, predecir si una casa es barata o cara. Con los datos disponibles, existen dos formas de comparar los modelos. Si bien las dos no tienen por qué dar los mismos resultados, son complementarias a la hora de tomar una decisión final.

models_cross = data.frame(
"modelo"= c('GLM','knn','Decision_Tree','Random_Forest','SVM'),
 "F_Medida_train" = c(F_reglog_train,F_knn_train,F_dt_train,F_rf_train,F_svm_train))
plot(x = models_cross$modelo, y= models_cross$F_Medida_train,fill=models_cross$modelo)

ggplot(data=models_cross, aes(x=modelo, y=F_Medida_train, fill=modelo)) + 
    geom_bar(stat="identity", position="dodge")

9 Curva ROC

# REGRESIÓN LOGÍSTICA
predictions_glm <- predict(model_glm, newdata = datos_train_rl, type = "response")
pred_glm <- prediction(predictions_glm, datos_train_rl$price_categ1)
perf_glm <- performance(pred_glm,"tpr","fpr")

# KNN
# model_knn=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=17, prob = TRUE)
predictions_Knn <- knn(scale(datos_train_numeric), scale(datos_train_numeric), cl=datos_train_limpio$price_categ1, k=17, prob = TRUE)
pred_knn <- prediction(attr(predictions_Knn,"prob"), datos_train_limpio$price_categ1)
perf_knn <- performance(pred_knn,"tpr","fpr")

# ARBOL DE DECISION
predictions_tree <- predict(decisiontree_model, newdata = datos_decision_tree, type = "prob")
pred_tree = prediction(predictions_tree[,2], datos_decision_tree$price_categ1)
pref_tree = performance(pred_tree, "tpr", "fpr")

# RANDOM FOREST
predictions_rf <- predict(randomforest_model, new_data=datos_train_limpio, type = "prob")
pred_rf <- prediction(predictions_rf[,2],datos_train_limpio$price_categ1)
perf_rf <- performance(pred_rf,"tpr","fpr")

# SVM
predictions_svm <- predict(modelo_svm, newdata=datos_train_limpio, probability = TRUE)
prob_svm<-attr(predictions_svm,"probabilities")
pred_svm <- prediction(prob_svm[,2], datos_train_limpio$price_categ)
perf_svm <- performance(pred_svm,"tpr","fpr")


plot(perf_glm,col="blue")
plot(pref_tree, col="red",add = TRUE)
plot(perf_rf,col="green", add = TRUE)
plot(perf_svm, col="yellow",add = TRUE)
plot(perf_knn, col="orange",add = TRUE)
legend(x="right", legend=c("GLM","DT","RF","SVM","KNN"), fill=c("blue","red","green","yellow","orange"), cex=0.8)

10 Conclusiones

  • nos quedamos con la acuracy ya que es la única medida que tiene en cuenta los aciertos de las casas caras y de las baratas (la F_medida no tiene en cuenta los TN y la sensibilidad y especificidad solo tienen en cuenta que acierte las caras o las baratas pero no que acierte en ambas)